//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

// The original of this code was found from:
//   http://www.netdenizen.com/emagnet/offaxis/elliptic_c.htm
//   http://www.netdenizen.com/emagnet/offaxis/iloopoffaxis.htm
//
// This code has been heavily modified.
//
// I got rid most of error detection, we loop a fixed level of times
// XXX need to test and investigate if this really really is ok

/*
 *
 * ELLIPTIC.C:  Functions for computing complete and incomplete elliptic
 *              integrals of the first and second kind, and an example
 *              of usage for computing magnetic fields due to a circular
 *              current filament.
 *
 * DISCLAIMER
 *
 * Although every effort has been expended to ensure correctness, this 
 * software is not guaranteed to be correct.
 *
 * You may do anything you want with this software, as long as you maintain
 * the existing comments intact, including the names of the original
 * authors.  If you modify this software, please indicate that fact in 
 * your comments.
 *
 * E. Dennison (Nov 17, 1998)
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>

#include "math3d.hpp"

/*
 * F(k): complete elliptic integral of the first kind
 */

#define DRF_C1 (1.0/24.0)
#define DRF_C2 (3.0/44.0)
#define DRF_C3 (1.0/14.0)
#define ITER 5

prec_t ellipse1(prec_t kk)
{
  prec_t xn = 0.0;
  prec_t yn = 1.0 - kk;
  prec_t zn = 1.0;

  for ( int i = 0; i < ITER; i++ )
    {
      prec_t xnroot = prec_t_sqrt (xn);
      prec_t ynroot = prec_t_sqrt (yn);
      prec_t znroot = prec_t_sqrt (zn);
      prec_t lambda = xnroot * (ynroot + znroot) + ynroot * znroot;
      xn = (xn + lambda) * 0.25;
      yn = (yn + lambda) * 0.25;
      zn = (zn + lambda) * 0.25;
    }

  prec_t mu = (xn + yn + zn) * (1.0 / 3.0);
  prec_t rmu = 1.0 / mu;
  prec_t xndev = 1.0 - xn * rmu;
  prec_t yndev = 1.0 - yn * rmu;
  prec_t zndev = 1.0 - zn * rmu;

  prec_t e1 = xndev * yndev;
  prec_t e2 = e1 - zndev * zndev;
  prec_t e3 = e1 * zndev;
  prec_t s = 1.0 + (DRF_C1 * e2 - 0.1 - DRF_C2 * e3) * e2 + DRF_C3 * e3;
  return s / prec_t_sqrt (mu);
}

/*
 * drd.c:   Compute the complete or incomplete elliptic integral of the
 *          second kind.
 *
 * Description:
 *
 *  For x and y non-negative, x+y and z positive, drf(x,y,z) = integral 
 *  from zero to infinity of 
 *
 *
 *            -1/2     -1/2     -3/2
 *  (3/2)(t+x)    (t+y)    (t+z)    dt.
 *
 *  If x or y is zero, the integral is complete.
 *
 *  *piErr returns non-zero if any arguments are invalid.
 *
 * Credits:
 *
 *  This function is adapted by E. Dennison from FORTRAN code written by:
 *
 *  Carlson, B.C.
 *  Notis, E.M
 *  Pexton, R.L.
 *
 */

#define DRD_C1 (3.0/14.0)
#define DRD_C2 (1.0/6.0)
#define DRD_C3 (9.0/22.0)
#define DRD_C4 (3.0/26.0)

prec_t ellipse2(prec_t kk, prec_t e1v)
{
  prec_t xn = 0;
  prec_t yn = 1.0 - kk;
  prec_t zn = 1.0;

  prec_t sigma = 0.0;
  prec_t power4 = 1.0;

  for ( int i = 0; i < ITER; i++ )
    {
      prec_t xnroot = prec_t_sqrt (xn);
      prec_t ynroot = prec_t_sqrt (yn);
      prec_t znroot = prec_t_sqrt (zn);
      prec_t lambda = xnroot * (ynroot + znroot) + ynroot * znroot;

      sigma = sigma + power4 / (znroot * (zn + lambda));
      power4 = power4 * 0.25;
      xn = (xn + lambda) * 0.25;
      yn = (yn + lambda) * 0.25;
      zn = (zn + lambda) * 0.25;
    }

  prec_t mu = (xn + yn + 3.0 * zn) * 0.2;
  prec_t rmu = 1.0 / mu;
  prec_t xndev = 1 - xn * rmu;
  prec_t yndev = 1 - yn * rmu;
  prec_t zndev = 1 - zn * rmu;

  prec_t ea = xndev * yndev;
  prec_t eb = zndev * zndev;
  prec_t ec = ea - eb;
  prec_t ed = ea - 6.0 * eb;
  prec_t ef = ed + ec + ec;
  prec_t s1 = ed * (-DRD_C1 + 0.25 * DRD_C3 * ed - 1.5 * DRD_C4 * zndev * ef);
  prec_t s2 = zndev * (DRD_C2 * ef + zndev * (-DRD_C3 * ec + zndev * DRD_C4 * ea));
  return e1v - kk * (sigma + power4 * (1.0 + s1 + s2) / (3 * mu * prec_t_sqrt (mu)));
}

prec_t ellipse2(prec_t kk)
{
  return ellipse2(kk, ellipse1(kk));
}

